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JING  Sfflt,  CHANGQING  HU*,  AND  CHI-WANG  SHU§ 


Abstract.  High  order  accurate  weighted  essentially  non-oscillatory  (WENO)  schemes  have  recently  been 
developed  for  finite  difference  and  finite  volume  methods  both  in  structured  and  in  unstructured  meshes. 
A  key  idea  in  WENO  scheme  is  a  linear  combination  of  lower  order  fluxes  or  reconstructions  to  obtain  a 
higher  order  approximation.  The  combination  coefficients,  also  called  linear  weights,  are  determined  by 
local  geometry  of  the  mesh  and  order  of  accuracy  and  may  become  negative.  WENO  procedures  cannot 
be  applied  directly  to  obtain  a  stable  scheme  if  negative  linear  weights  are  present.  Previous  strategy  for 
handling  this  difficulty  is  by  either  regrouping  of  stencils  or  reducing  the  order  of  accuracy  to  get  rid  of  the 
negative  linear  weights.  In  this  paper  we  present  a  simple  and  effective  technique  for  handling  negative  linear 
weights  without  a  need  to  get  rid  of  them.  Test  cases  are  shown  to  illustrate  the  stability  and  accuracy  of 
this  approach. 

Key  words,  weighted  essentially  non-oscillatory,  negative  weights,  stability,  high  order  accuracy,  shock 
calculation 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  High  order  accurate  weighted  essentially  non-oscillatory  (WENO)  schemes  have  re¬ 
cently  been  developed  to  solve  a  hyperbolic  conservation  law 

ut  +  V  ■/(«)=  0.  (1.1) 

The  first  WENO  scheme  was  constructed  in  [18]  for  a  third  order  finite  volume  version  in  one  space  dimension. 
In  [10],  third  and  fifth  order  finite  difference  WENO  schemes  in  multi  space  dimensions  are  constructed,  with 
a  general  framework  for  the  design  of  the  smoothness  indicators  and  nonlinear  weights.  Later,  second,  third 
and  fourth  order  finite  volume  WENO  schemes  for  2D  general  triangulation  have  been  developed  in  [4]  and 
[8].  Very  high  order  finite  difference  WENO  schemes  (for  orders  between  7  and  13)  have  been  developed  in 
[1].  Central  WENO  schemes  have  been  developed  in  [12],  [13]  and  [14]. 

WENO  schemes  are  designed  based  on  the  successful  ENO  schemes  in  [7,  23,  24],  Both  ENO  and 
WENO  use  the  idea  of  adaptive  stencils  in  the  reconstruction  procedure  based  on  the  local  smoothness 
of  the  numerical  solution  to  automatically  achieve  high  order  accuracy  and  non-oscillatory  property  near 
discontinuities.  ENO  uses  just  one  (optimal  in  some  sense)  out  of  many  candidate  stencils  when  doing  the 
reconstruction;  while  WENO  uses  a  convex  combination  of  all  the  candidate  stencils,  each  being  assigned 
a  nonlinear  weight  which  depends  on  the  local  smoothness  of  the  numerical  solution  based  on  that  stencil. 
WENO  improves  upon  ENO  in  robustness,  better  smoothness  of  fluxes,  better  steady  state  convergence, 
better  provable  convergence  properties,  and  more  efficiency.  For  a  detailed  review  of  ENO  and  WENO 
schemes,  we  refer  to  the  lecture  notes  [21,  22], 
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WENO  schemes  have  already  been  widely  used  in  applications.  Some  of  the  examples  include  dynamical 
response  of  a  stellar  atmosphere  to  pressure  perturbations  [3];  shock  vortex  interactions  and  other  gas 
dynamics  problems  [5],  [6];  incompressible  flow  problems  [26];  Hamilton- Jacobi  equations  [9];  magneto- 
hydrodynamics  [11];  underwater  blast-wave  focusing  [15];  the  composite  schemes  and  shallow  water  equations 
[16],  [17],  real  gas  computations  [19],  wave  propagation  using  Fey’s  method  of  transport  [20];  etc. 

A  key  idea  in  WENO  schemes  is  a  linear  combination  of  lower  order  fluxes  or  reconstructions  to  obtain 
a  higher  order  approximation.  The  combination  coefficients,  also  called  linear  weights,  are  determined  by 
local  geometry  of  the  mesh  and  order  of  accuracy  and  may  become  negative.  WENO  procedures  cannot 
be  applied  directly  to  obtain  a  stable  scheme  if  negative  linear  weights  are  present.  Previous  strategy  for 
handling  this  difficulty  is  by  either  regrouping  of  stencils  (e.g.  in  [8])  or  reducing  the  order  of  accuracy  (e.g. 
in  [12])  to  get  rid  of  the  negative  linear  weights.  In  this  paper  we  present  a  simple  and  effective  technique 
for  handling  negative  linear  weights  without  a  need  to  get  rid  of  them.  Test  cases  will  be  shown  to  illustrate 
the  stability  and  accuracy  of  this  approach. 

We  first  summarize  the  general  WENO  reconstruction  procedure,  consisting  of  the  following  steps.  We 
assume  we  have  a  given  cell  A  (which  could  be  an  interval  in  ID,  a  rectangle  in  a  2D  tensor  product  mesh, 
or  a  triangle  in  a  2D  unstructured  mesh)  and  a  fixed  point  xG  within  or  on  one  edge  of  the  cell. 

1.  We  identify  several  stencils  Sj,  j  =  1  such  that  A  belongs  to  each  stencil.  We  denote  by 

Q 

T  =  [J  Sj  the  larger  stencil  which  contains  all  the  cells  from  the  q  stencils. 

j  i 

2.  We  have  a  (relatively)  lower  order  reconstruction  or  interpolation  function  (usually  a  polynomial), 
denoted  by  pj(x),  associated  with  each  of  the  stencils  Sj,  for  j  =  1, ...,  q.  We  also  have  a  (relatively) 
higher  order  reconstruction  or  interpolation  function  (again  usually  a  polynomial),  denoted  by  Q(x), 
associated  with  the  larger  stencil  T. 

3.  We  find  the  combination  coefficients,  also  called  linear  weights,  denoted  by  71,  ...  ,  jq,  such  that 

Q(xG)  =  (1-2) 

3= 1 

for  all  possible  given  data  in  the  stencils.  These  linear  weights  depend  on  the  mesh  geometry,  the 
point  xG ,  and  the  specific  reconstruction  or  interpolation  requirements,  but  noton  the  given  solution 
data  in  the  stencils. 

4.  We  compute  the  smoothness  indicator,  denoted  by  0j,  for  each  stencil  Sj,  which  measures  how 
smooth  the  function  pj(x)  is  in  the  target  cell  A.  The  smaller  this  smoothness  indicator  pj,  the 
smoother  the  function  pj(x)  is  in  the  target  cell.  In  all  of  the  current  WENO  schemes  we  are  using 
the  following  smoothness  indicator: 

Pj=  E  /  |A| \a\-1{Dapj{x)fdx  (1.3) 

i<M<£ ' A 

for  j  =  1 ,  where  k  is  the  degree  of  the  polynomial  pj(x)  and  |  A]  is  the  area  of  the  cell  A  in  2D. 
This  factor  is  different  for  ID  or  3D:  the  purpose  of  it  is  to  bring  the  smoothness  indicator  invariant 
under  spatial  scaling. 

5.  We  compute  the  nonlinear  weights  based  on  the  smoothness  indicators: 


7 3 


(e  +  P j)'2 


(1.4) 
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where  7 j  are  the  linear  weights  determined  in  step  3  above,  and  e  is  a  small  number  to  avoid  the 
denominator  to  become  0.  We  are  using  e  =  10-6  in  all  the  computations  in  this  paper.  The  final 
WENO  approximation  or  reconstruction  is  then  given  by 

R(xG)  =  Y,UfPj(xG).  (1.5) 

i=i 

We  remark  that  all  the  coefficients  in  the  above  steps  which  depend  on  the  mesh  but  not  on  the  data  of 
the  numerical  solution,  should  be  computed  and  stored  at  the  beginning  of  the  code  after  the  generation  of 
the  mesh  but  before  the  time  evolution  starts. 

We  now  use  a  simple  example  to  illustrate  the  steps  outlined  above.  We  assume  we  are  given  a  uniform 
mesh  Ij  =  (./•,,  1 .  1  -j )  and  cell  averages  of  a  function  u(x)  in  these  cells,  denoted  by  u,;.  We  would 
like  to  find  a  fifth  order  WENO  reconstruction  to  the  point  value  u(xi+1/ 2),  based  on  a  stencil  of  five  cells 
{Ii-2,Ii-i , Ii,Ii+i , Ii+2},  with  the  target  cell  containing  the  point  xi+i/2  chosen  as  A  =  f,. 

In  step  1  above  we  could  have  the  following  three  stencils: 

Si  =  {!)_ a,/*— i,/j},  S2  =  {Ij-i,Ii,Ij+ 1},  S3  =  h+i ,  Ii+2}, 

which  make  up  a  larger  stencil 


T  —  {Ii-2,Ii-l,Ii,Ii+l,Ii+2}- 


In  step  2  above  we  would  have  three  polynomials  pj(x)  of  degree  at  most  two,  with  their  cell  averages 
agreeing  with  that  of  the  function  u  in  the  three  cells  in  each  stencil  Sj.  The  higher  order  function  Q(x )  is  a 
polynomial  of  degree  at  most  four,  with  its  cell  averages  agreeing  with  that  of  the  function  u  in  the  five  cells 
in  the  larger  stencil  T.  The  three  lower  order  approximations  to  u(xi+ i/2),  associated  with  pj(x),  in  terms 
of  the  given  cell  averages  of  u,  are  given  by: 


,  .  1_  7  _  11_ 

Pl(xi+l/'2)  ~  gU'-2  —  QUi~  1  +  ~QUi' 

,  ,  1_  5  _  1_ 

P2\Xi+ 1/2)  —  —QUi~  1  +  QUi  +  gW'j+H 

,  .  1  5  1 

P3\xi+l/2)  ~  7 A  +  QU'  +  1  ~  RU,:+2' 


(1.6) 
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Each  of  them  is  a  third  order  approximation  to  u(x $+1/2)-  The  higher  order  approximation  to  u(xj. 1-1/2)) 
associated  with  Q(x),  is  given  by: 


,  1  13  47  9  1 

Q\xi+ 1/2)  -  -  MUi~  1  +  Q{)u'  +  or\Ui+1  -  o^^Ui+'2, 


30  ‘  60 

which  is  a  fifth  order  approximation  to  u(xi+ 1/2)- 
In  step  3  above  we  would  have 


20 


20 


(1.7) 


1 


7i  = 


72  = 


10  5 

It  can  be  readily  verified,  using  (1.6)  and  (1.7),  that 


73  = 


10 


QGA+1/2)  =7i  /'1  1/2)  +72  P2(a:*+i/2)  +73  P3(xi+i/-2) 

for  all  possible  given  data  uj,  j  =  i  —  2,  i  —  1,  i,  i  +  1,  i  +  2. 
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Fig.  1.1.  Reconstructions  to  «(*!:+1/2)-  Solid  lines:  exact  function:  symbols:  numerical  approximations.  Left:  fifth  order 
WENO.  Right:  fifth  order  traditional. 


In  step  4  above  we  could  easily  work  out  from  (1.3)  the  three  smoothness  indicators  given  by 

13  _  1  o 

Pi  =  ( Ui-2  -  2u,_i  +  UiY  +  -  (Ui- o  -  4m-i  +  3 Wj): , 

13  1  ., 

/?2  =  Y2  (u*-i  “  2u*  +  W'*+l)“  +  ^  (Ui- 1  -  «>).]  )"  , 

13  2  1  _2 

P 3  =  12  Ui  ~  ^Ui+1  "*■  ^  (3«*  —  4uj+i  +  Uj+ o)  . 

We  notice  in  particular  that  the  linear  weights  71 , 72 , 73  in  step  3  above  are  all  positive.  In  such  cases,  the 
WENO  reconstruction  procedure  outlined  above  and  the  scheme  based  on  it  work  very  well.  In  Fig.  1.1  we 
plot  the  approximation  to  u(x )  for  a  discontinuous  function  u(x)  =  2x  for  x  <  0  and  u(x)  =  —20  otherwise, 
by  the  fifth  order  WENO  reconstruction  on  the  left  and  by  the  fifth  order  traditional  reconstruction  (1.7) 
on  the  right,  with  a  mesh  23  =  ( i  —  0.4965)A.t  with  A;c  =  0.02.  We  can  clearly  see  that  WENO  avoids  the 
over  and  undershoots  near  the  discontinuity. 

We  now  look  at  another  simple  example  where  some  of  the  linear  weights  in  step  3  above  would  become 
negative.  We  have  exactly  the  same  setting  as  above  except  now  we  seek  the  reconstruction  not  at  the  cell 
boundary  but  at  the  cell  center  23.  This  is  needed  by  the  central  schemes  with  staggered  grids  [12],  Thus, 
step  1  would  stay  the  same  as  above;  step  2  would  produce 

,  ,  1  1  23 

Pl(Xi)  =  —11;  2  + 

,  ,  1  _  13  _  1  _  „ 

P2(Xi)  =  i  +  j^Ui  ~  24«;+ 1-  l1-8) 

,  ,  23  1  1 

mxi)  =  +  i2Ui+1  - 

Each  of  them  is  a  third  order  reconstruction  to  u( 23).  The  higher  order  reconstruction  to  u( xp,  associated 
with  Q(x),  is  given  by: 

,  3  29  _  1067 _  29  _  3 

Cj(xi)  =  - Ui— 2 - Ui- 1  4 - Ui - u.i  +1  4 - UiM'i,  (1.9 

’  640  '  480  960  480  +  640  +  ; 
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Fig.  1.2.  Reconstructions  to  u(xi).  Solid  lines:  exact  function;  symbols:  numerical  approximations.  Left:  fifth  order 
WENO.  Right:  fifth  order  traditional. 


which  is  a  fifth  order  reconstruction  to  u(xj).  Step  3  would  produce  the  following  weights: 

9  49  9 

71  “  _ 80’  72  “  40’  73  “  _ 80’ 

Notice  that  two  of  them  are  negative.  The  smoothness  indicators  in  step  4  will  remain  the  same.  This 
time,  the  WENO  approximation,  shown  at  the  left  of  Fig.  1.2,  is  less  satisfactory  (in  fact,  even  worse  than 
a  traditional  fifth  order  reconstruction  show  on  the  right),  because  of  the  negative  linear  weights. 

We  remark  that  negative  linear  weights  do  not  appear  in  finite  difference  WENO  schemes  in  any  spatial 
dimensions  for  conservation  laws  for  any  order  of  accuracy  [10],  [1],  and  they  do  not  appear  in  one  dimensional 
as  well  as  some  multi-dimensional  finite  volume  WENO  schemes  for  conservation  laws.  Unfortunately,  they 
do  appear  in  some  other  cases,  such  as  the  central  WENO  schemes  using  staggered  meshes  we  have  seen 
above,  high  order  finite  volume  schemes  for  two  dimensions  described  in  [8]  and  in  this  paper,  and  finite 
difference  WENO  approximations  for  second  derivatives. 

While  on  approximation  alone  the  appearance  of  negative  linear  weights  might  be  annoying  but  perhaps 
not  fatal  (Fig.  1.2),  in  solving  a  PDE  the  result  might  be  more  serious.  As  an  example,  in  Fig.  1.3  we  show 
the  results  of  using  a  fourth  order  finite  volume  WENO  scheme  [8]  on  a  non-uniform  triangular  mesh  shown 
at  the  left,  which  has  negative  linear  weights,  for  solving  the  two  dimensional  Burgers  equation: 


(1.10) 


in  the  domain  [—2,2]  x  [—2,2]  with  an  initial  condition  Uo(x,y )  =  0.3  +  0.7 sin  (f(ir  +  y))  and  periodic 
boundary  conditions.  We  can  see  that  serious  oscillation  appears  in  the  numerical  solution  once  the  shock 
has  developed.  The  oscillation  eventually  leads  to  instability  and  blowing  up  of  the  numerical  solution  for 
this  example. 

The  main  purpose  of  this  paper  is  to  develop  a  simple  and  effective  technique  for  handling  negative  linear 
weights  without  a  need  to  get  rid  of  them.  Test  cases  will  be  shown  to  illustrate  the  stability  and  accuracy 
of  this  approach. 
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Fig.  1.3.  2D  Burgers'  equation.  Left:  non-uniform  triangular  mesh  used  in  the  computation.  Right:  fourth  order  WENO 
result  at  t  =  0.473,  cfl=0.2,  without  any  special  treatment  for  the  negative  linear  weights. 


2.  A  splitting  technique.  We  now  introduce  a  splitting  technique  to  treat  the  negative  weights.  It 
is  very  simple,  involves  little  additional  cost,  yet  is  quite  effective.  The  WENO  procedure  outlined  in  the 
previous  section  is  only  modified  in  step  5  in  the  following  way: 

5’  If  min(ji,  ...,7g)  >  0  proceed  as  before.  Otherwise,  we  split  the  linear  weights  into  two  parts: 
positive  and  negative.  Define 

it  =  \  tfi  +  0  |'y* |) ,  7“  =  7+  -  7 i,  *  =  h  ■■■,  Q  (2-1) 

where  we  take  6  =  3  all  the  numerical  tests.  We  then  scale  them  by 

7  t=lflo±,  i=l,...,q.  (2.2) 

j= i 

We  now  have  two  split  polynomials 

QH*G)  =  J27fPj(xG)  (2.3) 

3=1 

which  satisfy 


Q(xG )  =  a+Q+(xG)  -  o~Q-{xg).  (2.4) 

We  can  then  define  the  nonlinear  weights  (1.4)  for  the  positive  and  negative  groups  jf  separately, 
denoted  by  tof,  based  on  the  same  smoothness  indicator  Qj.  We  will  then  define  the  WENO 
approximation  /?±(a;G)  separately  by  (1.5),  using  and  form  the  final  WENO  approximation  by 

R(xg)  =  a+  R+(xg)  -  a ~  R~(xg). 


We  remark  that  the  key  idea  of  this  decomposition  is  to  make  sure  that  every  stencil  has  a  significant 
representation  in  both  the  positive  and  the  negative  weight  groups.  Within  each  group,  the  WENO  idea  of 


6 


Fig.  2.1.  WENO  approximations  with  the  splitting  treatment  for  negative  linear  weights.  Left:  approximation  to  u(xi). 
Right:  Burgers  equation,  solution  at  t  =  5/tt2  ,  cfl=0.2. 


redistributing  the  weights  subject  to  a  fixed  sum  according  to  the  smoothness  of  the  approximation  is  still 
followed  as  before. 

For  the  simple  example  of  fifth  order  WENO  reconstruction  to  u(xj),  the  split  linear  weights  correspond¬ 
ing  to  (2.1)  are,  before  the  scaling, 


,  9  __  9 

^  “  80’  ”fl  ~  40’ 


72 


49 

40’ 


9 

40' 


We  notice  that,  as  the  most  expensive  part  of  the  WENO  procedure,  namely  the  computation  of  the 
smoothness  indicators  (1.3),  has  not  changed,  the  extra  cost  of  this  positive/negative  weight  splitting  is  very 
small. 

However  this  simple  and  inexpensive  change  makes  a  big  difference  to  the  computations.  In  Fig.  2.1  we 
show  the  result  of  the  two  previous  unsatisfactory  cases,  the  fifth  order  WENO  reconstruction  to  u(xi)  in 
Fig.  1.2  left,  and  the  approximation  to  the  Burgers  equation  in  Fig.  1.3  right,  now  using  WENO  schemes 
with  this  splitting  treatment.  We  can  see  clearly  that  the  results  are  now  as  good  as  one  would  get  from 
WENO  schemes  having  only  positive  linear  weights. 

It  is  easy  to  prove  that  the  splitting  maintains  the  accuracy  of  the  approximation  in  smooth  regions. 
We  will  demonstrate  this  fact  in  the  following  sections.  We  will  also  demonstrate  the  effectiveness  of  this 
simple  splitting  technique  through  a  few  selected  numerical  examples  in  the  next  sections.  The  main  WENO 
schemes  we  will  consider  are  fifth  order  finite  volume  WENO  schemes  on  Cartesian  meshes,  and  the  third 
and  fourth  order  finite  volume  WENO  schemes  on  triangular  meshes.  In  both  cases  negative  linear  weights 
appear  regularly. 

The  calculations  are  performed  on  SUN  Ultra  workstations  and  also  on  the  IBM  SP  parallel  computer 
at  TCASCV  of  Brown  University.  The  parallel  efficiency  of  the  method  is  excellent  (more  than  90%). 


3.  2D  finite  volume  WENO  schemes  on  Cartesian  meshes. 

3.1.  The  schemes.  We  describe  two  different  ways  to  construct  fifth  order  finite  volume  WENO 
schemes  on  Cartesian  meshes.  Comparing  with  finite  difference  WENO  methods  [10],  finite  volume  meth- 
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ods  have  the  advantage  of  an  applicability  of  using  arbitrary  non-uniform  meshes,  at  the  price  of  increased 
computational  cost  [2]. 

We  define  the  cell: 


h,i  =  \xi-i,xi+i\  x  [yj_i,yj+i\  (3.1) 

for  i  =  1,  =  1,  ...,n,  Where  Ij  j  needs  not  be  uniform  or  smooth  varying. 

The  three-point  Gaussian  quadrature  rule  is  used  at  each  cell  edge  when  evaluating  the  numerical  flux 
in  order  to  maintain  fifth  order  accuracy.  Let  (xG ,yG)  denote  one  of  the  Gaussian  quadrature  points  at  the 
cell  boundary  of  Ijj  given  by  T  =  {.t  =  xi_i,yj_i  <  y  <  Vj+i}-  There  are  two  ways  to  perform  a  WENO 
reconstruction  at  the  point  (xG  ,  yG). 

Genuine  2D:.  The  first  WENO  reconstruction  is  genuine  2D  finite  volume.  We  can  see  that  there  are 
totally  nine  stencils  Sg.t  ( s,t  =  —1,0,1).  Each  stencil  Ssj  contains  3x3  cells  centered  around  /,:+sj+t. 
On  each  stencil  we  can  construct  a  Q 2  polynomial  (tensor  product  of  second  order  polynomials  in  x  and  y) 
satisfying  the  cell  average  condition  (i.e.  its  cell  average  in  each  cell  inside  the  stencil  equals  to  the  given 

i 

value).  Let  T  =  (J  SSt t,  which  contains  5x5  cells  centered  around  7,;j.  On  T  we  can  construct  a  Q 4 

S,t=  —  1 

polynomial  satisfying  the  cell  average  condition.  The  WENO  reconstruction  is  then  performed  according  to 
the  steps  outlined  in  sections  1  and  2. 

We  would  like  to  make  the  following  remarks: 

1.  By  using  a  Lagrange  interpolation  basis,  we  can  easily  find  the  unique  linear  weights. 

2.  Even  for  a  uniform  mesh,  a  negative  linear  weight  appears  for  the  middle  Gaussian  point  (xG  ,  yG)  = 
(Xj-i  ,yj)-  Such  appearance  of  negative  linear  weights  has  also  been  observed  in  the  central  WENO 
schemes  [12],  see  the  example  in  sections  1  and  2  before. 

3.  By  Taylor  expansions,  we  can  prove  that  the  smoothness  indicators  yield  a  uniform  fifth  order 
accuracy  in  smooth  regions.  See  [10]  for  the  method  of  proof. 

Dimension  by  Dimension:.  The  second  WENO  reconstruction  exploits  the  tensor  product  nature  of  the 
interpolation  we  use.  This  WENO  procedure  is  performed  on  a  dimension  by  dimension  fashion.  The  WENO 
schemes  applied  in  [5],  [6]  belong  to  this  class.  Consider  the  point  (x G,yG)  as  above.  First  we  perform  a 
one  dimensional  WENO  reconstruction  in  the  y  direction,  in  order  to  get  the  one  dimensional  cell  averages 
(in  the  x  direction)  w(»,yG).  Then  we  perform  another  one  dimensional  WENO  reconstruction  to  w  in  the 
x  direction,  to  obtain  the  final  reconstructed  point  value  at  (xG,yG). 

We  would  like  to  make  the  following  remarks: 

1.  For  a  scalar  equation,  the  underlying  linear  reconstructions  of  the  above  two  versions  are  equivalent. 
For  nonlinear  WENO  reconstructions  they  are  not  equivalent.  Both  of  them  are  fifth  order  accurate 
but  the  actual  errors  on  the  same  mesh  may  be  different,  see  Table  3.1  below. 

2.  For  systems  of  conservation  laws  such  as  the  Euler  equations  of  gas  dynamics,  both  versions  of  the 
WENO  reconstruction  should  be  performed  in  local  characteristic  fields. 

3.  The  dimension  by  dimension  version  of  the  WENO  reconstruction  is  less  expensive  and  requires 
smaller  memory  than  the  genuine  two  dimensional  version.  The  CPU  time  saving  is  about  a  factor 
of  4  for  the  Euler  equations  in  our  implementation.  The  computed  results  are  mostly  similar  from 
both  versions. 

In  the  following,  we  will  give  numerical  examples  computed  by  the  above  WENO  schemes.  Splitting 
technique  has  been  used  in  all  the  computations  when  negative  linear  weights  appear.  We  will  show  the 


Table  3.1 
2D  vortex  evolution 


genuine  FV 

dim-by-dim 

N 

A;r 

L°°  error 

order 

T°°  error 

order 

20 

6.71E-1 

4.38E-2 

5.26E-2 

40 

3.77E-1 

3.10E-3 

4.59 

5.66E-3 

3.86 

80 

2.01E-1 

1.20E-4 

5.15 

3.96E-4 

4.22 

160 

1.00E-1 

4.39E-6 

4.76 

7.96E-6 

5.62 

320 

5.00E-2 

1.88E-7 

4.53 

2.90E-7 

4.77 

results  for  both  smooth  and  discontinuous  problems. 

3.2.  2D  vortex  evolution.  First,  we  check  the  accuracy  of  the  WENO  schemes  constructed  above. 
The  two  dimensional  vortex  evolution  problem  [21],  [8]  is  used  as  a  test  problem. 

We  solve  the  Euler  equations  for  compressible  flow  in  2D 

Ut  +  f(U)x+g(U)y  =  0,  (3.2) 

where 

U  =  (p,pu,pv,E)rf 
f(U)  =  ( pu,pu 2  +  p,puv,u(E  +p))T, 
g(U)  =  (pv,puv,pv2  +  p,v(E  +  p))T . 

Here  p  is  the  density,  (u,  v)  is  the  velocity,  E  is  the  total  energy,  p  is  the  pressure,  related  to  the  total  energy 
by  E  =  7^-j-  +  kp{u2  +  v2)  with  7  =  1.4. 

The  setup  of  the  problem  is:  the  mean  flow  is  p  =  1,  p  =  1,  (u,v)  =  (1,1)  and  the  computational 
domain  is  [0, 10]  x  [0, 10].  We  add,  to  the  mean  flow,  an  isentropic  vortex  (perturbations  in  (u,v)  and  the 
temperature  T  =  £,  no  perturbation  in  the  entropy  S  =  ^-): 

(Su,Sv)  =  ^e°-5il~^(-y,x),  ST  =  -{l  £i-r^  6S  =  0. 

where  (x,y)  =  (x  —  5,  y  —  5),  r 2  =x2  +y2,  and  the  vortex  strength  e  =  5. 

We  use  non-uniform  meshes  which  are  obtained  by  an  independent  random  shifting  of  each  point  from  a 
uniform  mesh  in  each  direction  within  30%  of  the  mesh  sizes.  The  solution  is  computed  up  to  t  =  2.  Table  3.1 
shows  the  L°°  errors  of  p.  We  can  see  that  both  the  genuine  two  dimensional  finite  volume  WENO  scheme 
and  the  dimension  by  dimension  finite  volume  WENO  scheme  can  achieve  the  desired  order  of  accuracy 
while  the  genuine  two  dimensional  scheme  gives  smaller  errors  for  the  same  mesh. 

3.3.  Oblique  shock  tubes.  The  purpose  for  this  test  is  to  see  the  capability  of  the  rectangular  WENO 
schemes  in  resolving  waves  that  are  oblique  to  the  computational  meshes.  For  details  of  the  problem,  we  refer 
to  [10].  The  2D  Sod’s  shock  tube  problem  is  solved  where  the  initial  jump  makes  an  angle  0  against  the  x 
axis.  We  take  our  computational  domain  to  be  [0,  6]  x  [0, 1]  and  the  initial  jump  starting  at  (x,  y)  =  (2.25, 0) 
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Fig.  3.1.  Oblique  Sod's  problem.  Density  p.  Top:  contour,  genuine  two  dimensional  WENO;  middle:  contour,  dimension 
by  dimension  WENO;  bottom:  cut  at  the  bottom  of  the  computational  domain,  the  solid  line  is  the  exact  solution,  the  triangles 
are  the  genuine  two  dimensional  WENO  results,  and  the  circles  are  the  dimension  by  dimension  WENO  results. 


and  making  a  6  =  f  angle  with  the  x  axis.  The  solution  is  computed  up  to  t  =  1.2  on  a  96  x  16  uniform 
mesh.  In  Fig  3.1  we  plot  the  density  contours  computed  by  the  above  two  WENO  schemes  and  the  density 
cut  at  the  bottom  of  the  computational  domain.  We  can  see  that  both  schemes  perform  equally  well  in 
resolving  the  waves.  The  genuine  two  dimensional  scheme  gives  a  slightly  better  resolution  in  the  contact 
discontinuity  and  the  rarefaction  wave. 

3.4.  A  Mach  3  wind  tunnel  with  a  step.  This  model  problem  is  originally  from  [25].  The  setup 
of  the  problem  is:  The  wind  tunnel  is  1  length  unit  wide  and  3  length  units  long.  The  step  is  0.2  length 
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Fig.  3.2.  Forward  step  problem,  Ax  =  Ay  =  jg ,  ,  g|g  from  top  to  bottom.  30  contours  from  0.12  to  6.fl,  dimension 

by  dimension  WENO. 


units  high  and  is  located  0.6  length  units  from  the  left-hand  end  of  the  tunnel.  The  problem  is  initialized 
by  a  right-going  Mach  3  flow.  Reflective  boundary  conditions  are  applied  along  the  wall  of  the  tunnel  and 
inflow/outflow  boundary  conditions  are  applied  at  the  entrance/ exit.  The  corner  of  the  step  is  a  singular 
point  and  we  treat  it  the  same  way  as  in  [25],  which  is  based  on  the  assumption  of  a  nearly  steady  flow  in 
the  region  near  the  corner.  We  show  the  density  contours  at  time  t  =  4  in  Fig  3.2.  Only  the  results  from 
the  dimension  by  dimension  WENO  scheme  are  shown.  Uniform  meshes  of  Ax  =  Ay  =  Aj,  Aj,  _L-5  _L.  are 
used. 

3.5.  Double  Mach  reflection.  This  problem  is  also  originally  from  [25].  The  computational  domain 
for  this  problem  is  chosen  to  be  [0,4]  x  [0, 1].  The  reflecting  wall  lies  at  the  bottom,  starting  from  x  = 
Initially  a  right-moving  Mach  10  shock  is  positioned  at  x—  \,y  =  0  and  makes  a  60°  angle  with  the  x  axis. 
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Fig.  3.3.  Double  Mach  reflection.  Ax  =  Ay  =  (top  and  lower  left)  and  (middle  and  lower  right).  Genuine  two 
dimensional  WENO.  Blow-up  regions  at  the  bottom  for  details. 


For  the  bottom  boundary,  the  exact  post-shock  condition  is  imposed  for  the  part  from  x  =  0  to  x  =  ^  and 
a  reflective  boundary  condition  is  used  for  the  rest.  At  the  top  boundary,  the  flow  values  are  set  to  describe 
the  exact  motion  of  a  Mach  10  shock.  We  compute  the  solution  up  to  t  =  0.2.  Fig  3.3  and  Fig  3.4  show 
the  equally  spaced  30  density  contours  from  1.5  to  22.7  computed  by  the  genuine  two  dimensional  and  the 
dimension  by  dimension  WENO  schemes.  We  use  uniform  meshes  with  Ax  =  Ay  =  -jgg.  We  can  see 
that  the  results  from  both  schemes  are  comparable. 

4.  2D  finite  volume  WENO  schemes  on  triangular  meshes.  Both  third  and  fourth  order  finite 
volume  WENO  schemes  on  triangular  meshes  have  been  constructed  in  [8].  The  optional  linear  weights  in 
such  schemes  are  not  unique.  These  are  then  chosen  to  avoid  negative  weights  whenever  possible,  and  if  that 
fails,  a  grouping  (of  stencils)  technique  is  used  in  [8],  which  works  fairly  well  in  the  third  order  case  with 
quite  general  triangulation  but  can  yield  positive  weights  for  the  fourth  order  case  only  with  fairly  uniform 
triangulation.  In  this  section,  we  do  not  seek  positive  linear  weights  as  in  [8],  but  rather  use  the  splitting 
technique  to  treat  the  negative  linear  weights  when  they  appear.  For  scalar  equation,  the  scheme  is  stable 


12 


0.75 


Fig.  3.4.  Double  Mach  reflection,  Ax  =  Ay  =  ^jq  (top  and  lower  left)  and  (middle  and  lower  right).  Dimension  by 
dimension  WENO.  Blow-up  regions  at  the  bottom  for  details. 


in  all  runs.  For  systems  of  conservation  laws,  there  are  still  occasional  cases  of  overshoot  and  instability,  the 
reason  seems  to  be  related  to  characteristic  decompositions  and  is  still  being  investigated. 

4.1.  Accuracy  check  for  a  smooth  problem.  We  solve  the  2D  Burgers  equation  (1.10)  with  the 
same  initial  and  boundary  conditions  as  before  using  the  fourth  order  finite  volume  WENO  scheme  [8].  The 
solution  is  computed  up  to  t  =  when  no  shock  has  appeared.  The  meshes  used  are  1):  uniform  meshes 
with  equilateral  triangulation  and  2):  random  triangulation.  For  the  uniform  meshes  we  do  not  seek  positive 
weights  as  was  done  in  [8],  rather  we  use  the  splitting  technique  to  treat  the  negative  linear  weights  when 
they  appear.  Table  4.1  indicates  that  close  to  fourth  order  accuracy  can  be  achieved. 

4.2.  Discontinuous  problem  1:  Scalar  equation  in  2D.  Having  shown  the  stable  results  with  the 
splitting  treatment  of  negative  linear  weights  for  a  fourth  order  finite  volume  WENO  scheme  for  the  Burgers 
equation  in  section  2,  we  now  test  the  fourth  order  WENO  scheme  on  the  Buckley-Leverett  problem  whose 
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Table  4.1 

2D  Burgers  equation:  accuracy  check 


uniform  mesh 

nonuniform  mesh 

Ax 

L°°  error 

order 

Ax 

Z/°°  error 

order 

2.57E-1 

6.22E-4 

2.67E-1 

2.11E-3 

1.29E-1 

4.61E-5 

3.75 

1.26E-1 

2.35E-4 

2.92 

6.43E-2 

2.18E-6 

4.40 

6.32E-2 

2.90E-5 

3.03 

3.21E-2 

1.38E-7 

3.98 

3.34E-2 

2.61E-6 

3.78 

1.61E-2 

6.93E-9 

4.32 

1.66E-2 

2.71E-7 

3.24 

8.08E-3 

6.70E-10 

3.40 

7.44E-3 

1.57E-8 

3.55 

Fig.  4.1.  2D  Buckley- Lev erett  equation:  the  mesh. 


flux  is  non- convex: 


f(u) 


u2  +0.25(1  -u)2’ 


g(u)  =  0 


with  the  initial  data  u  =  1  when  —  4  <  x  <  0  and  u  =  0  elsewhere.  The  solution  is  computed  up  to 
t  =  0.4.  The  exact  solution  is  a  shock-rarefaction-contact  discontinuity  mixture.  The  mesh  we  use  here  is  a 
non-uniform  triangulation,  shown  in  Fig  4.1.  Fig  4.2  shows  that  the  waves  have  been  resolved  very  well. 


4.3.  Discontinuous  problem  2:  System  of  equations  in  2D.  We  consider  the  2D  Euler  equations 
in  the  domain  [—1, 1]  x  [0, 0.2].  The  Sod  and  Lax  shock  tube  initial  data  is  set  in  the  x  direction  and  periodic 
boundary  condition  is  applied  in  the  y  direction.  We  use  the  fourth  order  finite  volume  scheme  on  triangular 
meshes  to  solve  the  above  problem.  The  mesh  we  use  here  is  uniform.  But  we  do  not  seek  positive  weights  as 
was  done  in  [8],  rather  we  use  the  splitting  technique  in  section  2  to  treat  the  negative  linear  weights  when 
they  appear.  In  fact  we  set  deliberately  certain  linear  weights  to  be  negative  to  test  the  splitting  technique. 
Fig  4.3  shows  the  numerical  results  of  the  Sod  and  Lax  problems. 

It  seems  that  there  are  still  oscillations  and  instability  for  some  non-uniform  triangular  meshes  for  the 
fourth  order  WENO  schemes  applied  to  Euler  equations.  As  the  method  works  well  for  the  same  meshes 
with  a  scalar  equation,  the  problem  might  be  from  the  characteristic  decompositions.  This  is  still  under 
investigation. 


5.  Concluding  remarks.  We  have  devised  and  tested  a  simple  splitting  technique  to  treat  the  negative 
linear  weights  in  WENO  schemes.  This  technique  involves  very  little  additional  CPU  time  and  gives  good 
results  in  most  numerical  tests.  The  only  case  where  it  still  yields  oscillations  and  instability  is  when  a  fourth 
order  finite  volume  WENO  method  is  used  on  some  non-uniform  triangular  meshes  for  Euler  equations,  the 
reason  of  which,  presumably  related  to  characteristic  decompositions,  is  still  under  investigation. 
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Fig.  4.2.  2D  Buckley- Lev erett  equation  at  t=0.f,  with  splitting.  Left:  the  solution  surface;  Right:  the  cut  at  y  =  0.1  (solid 
line:  exact  solution,  symbols:  numerical  solution). 


Fig.  4.3.  Density  plot,  Left:  Sod  problem,  Right:  Lax  problem,  with  splitting.  Roughly  100  points  in  the  x  direction. 
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